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ABSTRACT 

The interpretation of gravity is valuable for understanding lithospheric plate 
motion and mantle convection. This paper compares postulated models of 
anomalous mass distributions in the earth with the observed geopotential as 
expressed in the spherical harmonic expansion. In particular, models of the 
anomalous density as a function of radius are found which can closely match 
the average magnitude of the spherical harmonic coefficients of a degree L 
These models include (1) a two-component model consisting of an anomalous 
layer at ~200 km depth (below the earth’s surface) and at ~1500 km depth 
(2) a two-component model where the upper component is distributed in the 
region between 1000 and 2800 km depth, (3) a model with density anomalies 
which continuously increase with depth more than an order of magnitude. 
Models of anomalous density spread into laterally extended caps or confined 
to depths above 1000 km do not match the properties of the observed field 
unless very special and unlikely arrangements of anomalous densities are 
assumed. If the average coefficient in a degree of the geopotential is a con- 
sequence of anomalous densities which vary as a function of radius, then 
only models in which the anomalous densities increase with depth and extend 
to great depth, probably to the core-mantle boundary, appear able to match 
the geopotential. The most reasonable geophysical interpretation of deep 
lateral density variations is that they are associated with lower mantle 
convection. 
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LATERAL DENSITY ANOMALIES AND THE EARTH’S 


GRAVITATIONAL FIELD 

I. INTRODUCTION 

An unresolved problem of geophysics is the motive force of plate tectonics. Since the motion 
of plates is strongly correlated with earthquakes, information which may help to understand 
this motion is useful. While it is generally agreed that plate motion is associated with con- 
vective motion, the question of whether this convection is confined to the upper mantle or 
arises from the lower mantle is controversial at present. It is the purpose of this paper to 
study the relation of mass anomalies to the geopotential in order that the results may be 
used to support, refute or aid in developing geophysical hypotheses concerning mantle 
motions. 

The lateral density anomalies are of significance in studying mass imbalances due to convec- 
tion. The radial distribution of density contributes only to the central force term of the geo- 
potential and is obtained from a combination of geophysical and astronomical data and 
Monte Carlo calculations (Haddon and Bullen, 1969; Press, 1970; Gilbert, Dziewonski and 
Brune, 1973). These computations provide a density p(r) to which anomalous densities are 
referred. These lateral density anomalies A p do not alter the amount of mass in a spherically 
symmetric shell when summed over the shell. 

Lateral density anomalies are expressed in the coefficients of the terms above the first degree 
and order in the spherical harmonic expansion of the geopotential. Satellite perturbation 
data has been used to derive an extensive set of coefficients for the spherical harmonic ex- 
pansion of the geopotential. These coefficients, often termed “C’s and S ’s” are particularly 
convenient to extract from long arc satellite orbits. Often, there are individual coefficients 
which are resonant with a particular satellite’s orbit and may be individually determined 
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(Wagner, 1976; King-Hele, Walker and Gooding, 1974, 1975; Klokocnik, 1975). More gen- 
erally, geopotential field models are constructed by simultaneously solving for the coefficients 
from large data bases constructed from many satellite orbits (Lerch et al„ 1977; Gaposchkin, 
1974). The satellite arcs include a variety of inclinations, periods and mean distances from 
the earth and therefore yield a data base which spans the globe. A global distribution is neces- 
sary to make a meaningful determination of the coefficients as they are computed from a 
global integration. The global knowledge allows the possibility of separating global contribu- 
tions from local contributions and of estimating the structure of lateral density variations as 
a function of depth. 

But although the spherical harmonic expansion contains global geophysical information, the 
extraction of this geophysical information is not a straightforward procedure. The major dif- 
ficulty is that knowledge of the the external potential does not allow a unique mathematical 
inversion procedure to recover internal density variations. An additional difficulty is that the 
density variations within the earth probably are not distributed in the shape of any individual 
harmonic with the exception of the hydrostatic bulge of the earth. The hydrostatic bulge is 
mainly expressed in the second (zonal) harmonic, constituting the overwhelming proportion 
of it, and 1/4 of the fourth zonal harmonic. The hydrostatic contribution is well understood 
(Jeffreys, 1963), and this paper considers possible models of the nonhydrostatic part of the 
earth’s potential using the method of calculating the nonhydrostatic spherical harmonics 
from assumed density distribution. 

There have been several previous investigations of the geophysical implications of the spheri- 
cal harmonic coefficients of the earth’s gravity field. Allan (1972), using a method of analy- 
sis initially described by Cook (1967) and employed by Guier and Newton (1965), and Khan 
(1970), concluded that the earth’s potential could be modeled in terms of two components 
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where the harmonics below degree 6 or 8 could be described by density 
anomalies at an effective depth below the earth’s surface of ~1700 km, 

and higher harmonics at a depth of ~300 km. Khan (1977) and Lambeck 
(1976), using essentially similar techniques of analysis concluded that all harmonics could 
be described by anomalies distributed between 0 and 800-1000 km depth. Suggested 
sources of the earth’s gravitational field also include the core-mantle interface (Hide and 
Horai, 1968; Hide and Malin, 1970; Khan and Woolard, 1970), the upper mantle transitions 
zones (Bott, 1971), and the upper mantle (Higbie and Stacey, 1970). Occasionally attempts 
have been made to interpret individual harmonic coefficients with particular features (Khan 
and O’Keefe, 1974); Wang, 1966). McDonald (1963) suggested that the non-hydrostatic 
portion of the second harmonic was a fossil bulge remaining from an earlier rate of rotation 
of the earth. 

This paper utilizes the approach employed by Goldreich and Toomre (1969) and by Pollack 
(1973): a mass distribution is postulated and the resulting coefficients of the spherical har- 
monics are calculated. Whereas Goldreich and Toomre calculated only the effect of random 
mass distributions on the second harmonics, this paper calculates the harmonics up to 
degree and order 25. This approach investigates the mathematical implications of converting 

a lateral anomalous density into spherical harmonics first and then compares the results with 
actual field models and other geophysical data or theory. In this way, no geophysical 
assumptions are required initially (except that the density anomalies are small in comparison 
with the density) and a number of plausible geophysical configurations are investigated. 
These include mass points, plumes, convection rolls, extended caps and randomly distributed 
masspoints covering layers on the globe. The variation in anomalies with depth is studied in 
the two limiting cases; assuming the anomalies to be in a narrow layer-single depth— and 
assuming the anomalies to be distributed throughout most of the earth’s mantle. An 
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important difference in this work is that the results of simulation are compared with observed 
average coefficient values in each degree, rather than basing conclusions on a least square fit 
which has been summed over all degrees. The ability of a model to satisfy all degrees equally 
provides a more stringent criterion to evaluate the model. 

Procedure 

The external potential of the gravity of the earth satisfies LaPlace’s Equation (MacMillan, 
1958) 

V 2 V = 0 (1) 

The external potential may be expressed as 

V = v 1 + v 2 

where Vj = GM/r is the potential of a central force. The anomalous external potential V 9 is 

°° g a £ 

V 2 = GM ^2^2 ~ P 5 „,0 in cos mX + s sm sin mX) (2) 

£=0 ni=0 * 

where G is the Newtonian gravitational constant, M the total mass of the earth, a e the radius 
of the earth, 0 the latitude, X the longitude, r the distance from the center of the earth and 
P em are normalized spherical harmonics or associated Legendre polynomials. The C Cm and 
S fini are estimated from data such as satellite perturbations, altimetry and ground-based data. 
The normalized coefficients are more suitable for geophysical interpretation as the unnormai- 
ized coefficients are related to the normalized coefficients by a factor of p2(2C + 1 ) ^ 

L (£ + m)!j 

1/ 

if m =£ 0 and (2£ + 1 ) n if m = 0 and thus the unnormalized tesserals decline rapidly in a given 
degree £. This effect in the unnormalized coefficient along with the circumstance that only 
the zonals are of any import in the hydrostatic field, has in the past sometimes caused the 
tesserals to be considered of less importance. 
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If a density distribution due to a set of mass points m. is assumed, the resulting spherical 
harmonic coefficients can be calculated readily: 


c 8m = Z m, (2j7T) Ptm(sin0i)COSmX ‘ 

S em = L m, ^ !) < sin ^ sinm \ 


( 3 ) 


where nij is a “mass point,” normalized to the earth’s mass and located at depth dj below 
the surface of the earth, latitude <f> i and longitude X., The point mass form is most conven- 
ient numerically but the mass may be spread in a volume element having density Apj in which 
case the above equation would take the form of an integral rather than a summation. The 
coefficients of the potential of a unit spherical cap on the earth’s surface are 



P c+ j (cos a.) - P (cos a) _ 

7 1 P Cm («n0) 

(2fi + l) 2 (1 - cos a) 


cos mX 


sin mX 



where a is the generating angle of the cap (Pollack, 1973). The computation is generalized 

/% - d 

for a mass m at a depth d by the factor m { 

V a e 

The amplitude of a degree is often used to make correlations with geophysical properties. 
The quantity e 

°l = 23 < C Sm + ^«m) (S > 

m=0 

is often termed the “degree variance” or the “power spectrum” or the “spectral energy.” 
This dimensionless quantity is rotationally invariant. It is related to the root mean 
square coefficient variation (C^ or Sg m > in a degree fi: 

°e ~ (2fi + 1) <C B 2 m .S c 2 m > (6) 
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Kaula (1966) observed that the harmonics follow this empirical rule: 

<C Sm ,S- Sm > - lQ-*/8» (7) 

Pollack (1973 noted that the degree variance of a unit mass point on the surface of the earth 
may be normalized by multiplying by (2 8 +1). This motivated his definition of the quantity 

*, = o s (2£+l) i4 (8) 

where a mass point of unit value of the surface of the earth will have <f> £ (8) = 1 at all degrees. 
This quantity is here termed the “relative amplitude” since it reflects the amplitude of a de- 
gree computed from an arbitrary mass distribution relative to the amplitude of a degree com- 
puted for a unit surface mass point. The relative amplitude is used in this paper to study the 
spectral properties of the harmonics because it is convenient for geophysical interpretation 
and graphical presentation. As Pollack showed, <i> e computed for a given total mass will 
diminish with degree when the mass is lowered from the surface or spread into a regular ex- 
tended feature. A field arising from non-organized crustal density distributions will exhibit 
a flat behavior of d? c as 8 increases. 
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fl. RESULTS OF SIMULATION OF STOKES COEFFICIENTS 
FROM DENSITY MODELS 

This section shows the results of actual numerical experiments obtained from postulating vari- 
ous mass distributions and calculating sets of coefficients in the spherical harmonic expansion 
of the geopotential (Eq. 3). The resulting coefficients or the general properties of a set of 
coefficients are compared with the actual coefficients measured in satellite geodesy. This is 
a fruitful approach because the associated Legendre polynomials in Eq. 3 are too complica- 
ted to allow an intuitive grasp of the effect of a mass configuration on the spherical har- 
monic coefficients and because it is possible to develop general statistical properties of a set 
of coefficients. It is not possible with this approach to infer a unique mass distribution but 
it is possible to evaluate the ability of interesting or geophysically plausible mass configura- 
tions to explain the gravitational field. The configurations used for simulations include a 
mass point, a vertical column of mass points, an extended spherical cap, sets of random mass 
points globally distributed at various depths, sets of concentric caps and convection rolls. 

Modern fields are now obtained to high degree and order. The relative amplitudes ( ( !> e from 
Eq. 8) for two field models, GEM-9 (Lerch et ah, 1977) and SAO-II (Gaposchkin and Lam- 
beck, 1 97 1) are shown in Fig. 1 and Fig. 2. It is noted that this quantity declines by a 
factor of 4 between degree 3 and degree 8 and declines less precipitously in higher degrees. 

a. Mass Point as a Function of Depth 

The spherical harmonic coefficients are simulated for a disturbing field composed of one 
mass point. This particularly simple distribution is able to demonstrate interesting general 
properties. The behavior of the relative amplitude as a function of degree is shown in Fig. 3 
for mass points at three depths, 0 km, 500 km and 1 500 km. Pollack (1973) showed that a 
point mass produces the ‘richest spectrum’ when the relative amplitude is compared with that 
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Figure 1 . Relative Ampltude of the GEM-9 Field Model 
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Figure 2. Relative Amplitude of the SAO II Field Model 


9 




<J> X 10 



Figure 3. Effect of Depth on Relative Amplitude 
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arising from extended sources such as a spherical cap or a spherical rectangular block. Also, 
if a density anomaly is less than 5° in lateral extent, the decline of <h fi is not substantially dif- 
ferent from to a point mass in the harmonics below degree and order 30. Therefore, al- 
though the mass point is not a completely realistic geophysical model, it is useful in describ- 
ing the contribution of a short wavelength density anomaly. 

If a mass point is placed on the surface of the earth, the relative amplitude is constant and 
equal to the value of the mass in earth mass (em) units. Thus a mass point of 8. X 10~ 6 em 
produces a value of <h e = 8. X 10 -6 em for all fi. If the mass point is placed slightly below the 
earth’s surface, then 4> c will decline with increasing degree and even the contribution to 
degree two will be diminished. At a depth of 500 km, a mass of 8. X 10 -6 em gives 
«I> 2 = 7. X 10 -6 , at 1500 km $ 2 = 5. X 10~ 6 . This effect increases sharply with increasing 

K - h V 

degree, according to the factor l ! in eq. 3 so that deep anomalies, as for example at 

\ 3 e / 

the core-mantle interface, will not contribute to more than the low degree harmonics, if at 
all perceptible in the earth’s measured field. 

In comparing the decline of the relative amplitude of the low degree terms with that of a real 
field, the decline of the real field up to degree 12 is followed most closely by a mass point at 
1500 km depth. After degree 13, the real field does not decline as steeply as a mass point at 
1000 km depth. But, below degree 12, the real field cannot be matched by a point mass 
located above 1000 km depth. This result derived from a mass point assumption suggesting 
that the field arises from a source located below 1000 km depth is consistent with earlier 
findings by Allan (1972), who used similar assumptions. The magnitude of a mass point 
matched to a low degree field originating at a depth of 1500 km is 1.7 X 10~ 5 em (1 em = 

6 X 10 24 kg) located at depth = 1500 km (Fig 4). If placed at 1000 km depth, this mass 
does not match except for the second harmonic. 
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Figure 4. Relative Amplitude From a Mass Point at 1 ,000 and 1 ,500 km Depth 

(m - 1.7 X 10" 5 em) 
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b. Mass Point as a Function of Geographic Location 

The properties of the relative amplitude discussed above are independent of the geographic lo- 
cation of the mass point (by the rotational invariance of spherical harmonics); however; the 
geographic location of a mass point on a sphere controls the distribution of the amplitude in a 
degree £ among the m orders. There is a particularly striking effect when the point mass is 
located near the poles or in a high latitude; the zonal and low order tesserals contain all the 
magnitude and the higher order tesserals and sectorial contribute negligibly to the degree. 
Table 1 contains a list of the C x 2 m computed for a mass anomaly located near the pole: 
the magnitude of the coefficients declines strongly with increasing order. Similar effects are 


Table 1 

Twelfth Degree Field from a Polar Mascon 


m 

C 12-m 


1 

0.22 X 10“ 06 


2 

-0.34 X 10~ 06 


3 

0.22 X 10“ 06 


4 

-0.93 X 10“ 07 

LAT = 76° 

5 

0.29 X 10" 07 

LONG = 180° 

6 

-0.71 X 10“ 08 

HEIGHT = 0. 

7 

0.13 X10" 08 

MASS = -y 4 X 10“ 6 em 

8 

-0.22 X 10“ 09 


9 

0.28 X 10“ 10 


10 

-0.29 X 10“ 11 


11 

0.22X10“ 12 


12 

-0.11 X 10“ 13 
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apparent at all degrees. Mass points near the equator produce no decline in the normalized 

coefficients of the potential field. Table 2 shows the values of C 12 m for a mass point located 
coeff 

near the equator. These equal values of the coefficients for all orders m is compatible with the 
behavior of the real field. The discrepancy between the behavior of a polar and equatorial mass 
distribution is a direct consequence of the numerical behavior of the normalized Legendre 
polynomials: Fig 5 show P 30 > P 3 1 » P * 32 * and ^33 as a function of latitude. This decline in 
amplitude as a function of order as the latitude increases becomes more pronounced as the 
degree increases. A field arising from only density anomalies located at high latitudes would 
be a mainly zonal field. 


Table 2 


Twelfth Degree Field from an Equatorial Mascon 


m 

^T 2 -m 


1 

-0.14 X 10" 07 


2 

0.49 X 10 ~ 08 


3 

0.74 X 10 " 08 


4 

-0.16 X 10 ~ 08 

LAT = 20° 

5 

0.22 X 10 “ 08 

LONG - 200° 

6 

-0.69 X 10 ” 08 

HEIGHT = 0. 

7 

-0.60 X 10 ~ 08 

MASS = % X 10 " 6 em 

8 

0.17 X 10 ' 07 


9 

-0.33X10" 08 


10 

-0.17 X 10 -07 


11 

0.18 X 1Q ~ 07 

.. . . 


12 

-0.67 XI 0 " 08 



14 






AMPLITUDE 


.6 


.4 


.2 


0 


-.2 


-.4 


•.6 

- 80 ° - 60 ° - 40 ° - 20 ° 0 ° 20 ° 40 ° 60 ° 80 ° 

LATITUDE 

Figure 5. Third Degree Legendre Polynomials vs. Latitude 
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These results should be considered in discussions of the source of “the third harmonic” 
(Khan and O’Keefe, 1974; Kaula, 1972); it is pointed out here that an “Antarctic anomaly” 
invoked to explain the high value of C 30 = .95 X 10” 6 cannot simultaneously explain the 
high value of S 33 - 1.4 X 10~ 6 . In fact, numerical explorations failed to locate any point 
on the globe which could satisfactorily represent all the third harmonics, including the 
largest, the C 31 = 2.0 X 10~ 6 . 

c. Column and Random Layer Distributions 

The effect of distributions extending through a range of depths was investigated using a local 
distribution arranged in a column of point masses and a global distribution consisting of 
mass points with random magnitudes inserted into equal area blocks. The local anomaly 
was studied using a column of masses at the same geographic location placed at depths of 
1 00 km apart. The global distribution was investigated by constructing layers separated by 
100 km in depth composed of random valued mass points placed in equal area blocks. 

The simplest construction of reflecting anomalies at many depths is a column of equal- 
valued mass points. A column of 1 1 mass points, located at the same geographic position 
(latitude — 20°, longitude = 200°) but varying in depth location by 100 km, was constructed 
placing a mass point of .5 X 10~ 6 em at 0 km, - 100 km, -200 . . - 1000 km. The relative 
amplitude resulting from this configuration is the smooth line shown in Fig. 6. 

Similar results were obtained from a set of 1 1 global layers placed at 0, - 100,-200, 

- 1000 km. Each infinitesimally thin layer was constructed by dividing the spherical surface 
into equal-area 10° X 10° blocks (Paul, 1973) and inserting a random mass point into each 
block. The random value of the mass point was generated from a Gaussian distribution 
taken about 0 with a standard deviation of 2 X 10“ 8 earth masses. The results are shown as 
squares in Fig. 6. These results are compatible with those from the column of mass points, 
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Figure 6. Relative Amplitude from Anomalies Above 1,000 km Depth 


17 




except for a slight offset of the curve due to a slightly different effective mass. The global 
layer of mass points, which is in essence modeling short wavelength uncorrelated anomalies, 

produces relative amplitude that scatters about the smooth line resulting from an individual 

■ , i ; 

column, In comparing with a field mddel, indicated in Fig. 6 by circles, the simulated 
energy spectra cannot reproduce the sharp decline in the low degrees. Even if the harmonics 
below four are deleted, the simulated amplitudes do not decline as rapidly as the observed. 
Further, it was found that random layers of equal anomalous densities which begin at the 
surface and extend to depths down to the core-mantle interface do not produce substantially 
different fields from those extending to only 1000 km. Thus an anomalous field of uncorre- 
lated short wavelength features with equal mass per spherical area at each depth would not 
usually produce the observed field. 

A single global layer of random mass points at a depth of 1500 km does replicate the decline 
of the relative amplitude of the real field in the low orders as shown in Fig. 7. The layer was 
constructed by placing 206 mass points averaging 1 X 10“ 6 em in 10° X 10° global equal 
area blocks. The overall decline of the $ is well-matched but these differences 
remain 1) the random field produces more scattered values of the low degree harmonics 
2) the observed high degree harmonics (2 > 12) do not vanish as do the modeled, a differ- 
ence corrected by adding an upper component similar to the ‘two-component’ model ob- 
tained by Allan (1972). 

The assumption of two narrow layers of density anomalies imposed severe geophysical con- 
straints, and so a modification of the two-component model was investigated, where the 
low degree component was spread in layers between 1000 and 2800 km depth. Fig. 8 shows 
the relative amplitude of a two-component model assuming masses of 1, X 10 -6 em arranged 
in a column of 19’points at 1000, 1100, 1200 ... 2800 km depth and one mass point at 
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Figure 7. Relative Amplitude from an Anomalous Layer at 1 ,500 km Depth 
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200 km depth. (The high degree component could just as well be modeled by spreading thp 
density anomaly at 200 km depth through several hundred km because the relative 
amplitude does not resolve between these two possibilities.) The $ of this modified 
tworcomponent model matches that of the observed field well. This model implies a greatly 

i .... ; 

increased anomalous density for the lower component because the deep mass points are 
placed every 100 km apart whereas there is only one mass point in the upper 1000 km. 

Because the modified two-component model places possibly artificial constraints on the 
density variations, a model employing monotically increasing density anomalies was tested. 

In this model, an effective mass of .25 X 10 -6 em was placed on the surface layer and a 
mass incremented by .025 X 10~ 6 em was placed on each successive layer deeper by 100 
kni, that is, m* = .250 X 10- 6 ‘at 0. km, .275 X 10" 6 at 100 km, .300 X 10" 6 at 200 km, 

.325 X 10" 6 at 300 km. . .... . .95 X 10 -6 at 2800 km depth. This produced relative ampli- 

. ( .. 

tudes (Fig. 9) with properties similar to that observed, and a reasonable match in magnitudes. 
Veiry close matches in magnitude can be obtained by any one of a number of possible 
adjustments in density variation. Also, discontinuities in anomalous density such as are 
known to occur in the radial distribution of density can reproduce the observed energy spec- 
trum since the spherical harmonic coefficients do not resolve this type of model differences. 
The mass points vary by a factor of 4 from the surface to the core-mantle interface, Consid- 
ering that the area of a layer shrinks by a factor of 4 in this distance, this model implies 
that the anomalous density increases by a factor of 16 at the bottom of the mantle. 

In order to investigate the origins of the second harmonic, the percentage contribution to 
the $2 (amplitude of the second degree harmonic) from each depth can be computed (Fig. 
10). Each of the 29 layers contributes between 2.5 and 3.5% of the total value of re- 
sulting in a nearly flat curve with a broad maximum between 1000 and 2000 km depth. 

21 







2000 km depth. This monotonically increasing density model provides no single layer or 
region which can be said to produce the second harmonics. 

The scatter produced by the field composed of individual mass points is a function of the 
number of anomalies involved. Investigations with a set of 827 mass points inside 5° X 5° 
equal area blocks and 92 mass points inside 15° X 15° blocks produced results in accord 
with those expected from considering the problem as a random walk process. The scatter 
is highest in the low degrees as a statistical consequence of summing over fewer values of 
C _ and S „ . As the number N of mass points increases, the net values of <I> 0 tend to in- 
crease as^/W. A set of N masses of nij each on a layer yield an effective mass equal to the 

Nnij 

total mass divided by>/N, that is, m eff =— ^= V /Nm.. Vertical inhomogenieties combine 
in the same manner, except that the strength of a layer is diminished according to the depth. 

Extended Sources 

The relative amplitude 3> e is not only a function of the depth of an anomalous density, it is 
also a function of geographic extension of the anomaly. Pollack (1973) showed the relative 
amplitude of a degree resulting from spherical caps or rectangles on the surface of the earth 
declines rapidly with degree if the anomalous source covers a significant area on the surface. 
The slope of the decline increases as the area of the lateral extension increases. The impact 
of a cap on is computed, where a cap is defined as a circular area on a shell generated by 
a central angle a. For a cap of radius a = 5°, <J>g is nearly like a point mass; after £ = 20, 
was reduced to about 70% of that of an equivalent mass point. However, from a 10° cap 
was reduced over an order of magnitude by £ = 20. As the radius of the cap increases, the 
decline of <I>g is no longer monotonic; a strong cycloidal behavior occurs, although the oscil- 
lations are still superimposed on a declining trend. The amplitude of a 30° cap has declined 
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an order of magnitude by £ = 7, and of a 60° cap by degree 3. Pollack also found that mass 
anomalies in the shape of a rectangle on a spherical shell produced similar results for 

The computations in this study confirmed Pollack’s results and considered spherical caps 
with more structure than uniform density. By constructing a set of spherical caps with a 
common center but with varying radii, a layered density can be simulated. These layered 
cap with mass anomalies increasing toward the center are used to model features observed 
in the real gravity field. 

Two models of layered caps were constructed for the Antarctic low using the geoidal 
anomalies shown in Khan and O’Keefe (1974). One cap was modelled after the anomaly 
referred to the reference ellipsoid (flattening f = 1/298.255) and one cap to an ellipsoid with 
hydrostatic flattening (f = 1/299.75). The model anomalous to the reference ellipsoid is 
centered at 0 = 72° S, X = 183°E. It is composed of a set of superposed concentric caps 
having m A = - 1.3 X 10“ 6 em in a 20° cap, -.6 X 10“ 6 em in a 14° cap, -8.3 X 10“ 6 ’em in a 
10° cap, - .2 X 10“ 6 in a 7° cap and - .1 X 10“ 6 'em in a 4° cap. The model referred to the 
ellipsoid with a hydrostatic flattening is centered at 0 = 76°S, X = 177°E, and consists of 
caps of size 26° with-1.3 X 10“ 6 em, 18.2° with-.6 X 10“ 6 em, 11.3° with -.3 X 10“ 6 
em, 8.8° with -.2 X 10" 6 em and 6.2° with-.l X 10“ 6 em. All caps were assumed to be on 
the surface of the earth. The 4> e calculated for both models dropped an order of magnitude 
in the first 20 degrees (Fig. 1 1). The two models produced similar curves, except that the 
hydrostatic flattening model varied slightly faster and showed the beginnings of a cycloidal 
shape in the first 25 degrees. An extended feature located at a depth would show a steeper 
decline in the relative amplitude. If an extended feature contains an intense central anomaly 
relative to the fringes, the relative amplitude will tend toward that of a mass point. 
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Thus, regular extended features near the earth’s surface demonstrate a sharp decrease in $ with 
degree. Further, this decrease differs from that of the real field and of fields simulated from 
depth in that the decline is convex rather than concave. Therefore, this type of feature does 
not readily match the observed field; further, short wavelength anomalies can be expected 
to predominate as the degree of the field increases. 

Organized surface elements may have unpredictable amplitudes. Convection rolls pro- 
vided a motivation for simulations of an organized type of extended source. The possibility 
of convection rolls underneath the Pacific plate was calculated from theoretical considera- 
tions by Richter and Parsons (1975) and tentatively identified in the observed field above 
degree 12 by Marsh and Marsh (1976). This type of anomaly was modeled with 4 longitudi- 
nal rows separated by 6° of longitude and spanning 9° of latitude. Each row consisted of 10 
mass points of magnitude ± 1 1 X 10~ 6 I em each. Alternate rows on the earth’s surface 

from 0 — 20 to 29° contain positive mass points at X = 150°, negative masses at X = 156°, 

. - m 

positive at^j = 162° and negative at X = 168°. The relative amplitude of this organized ele- 
ment is markedly different from random or extended features. The resulting shown in 
Fig. 12 increases non-monotonically with high degree. The expected average { I> 2 from a set 
of 40 random mass points of 1 X 10 -6 em on the surface is ~6. X 10" 6 for all degrees. The 
3> c for the modeled convection roll is 14. X 10“ 6 at degree 25, over twice that of an equiva- 
lent set of random mass points. Thus organized elements in the earth’s crust may contribute 
significantly to high degrees in a definitely non-intuitive manner. 
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DEGREE NUMBER 

Figure 1 2. Relative Amplitude of a Convection Roll 





III. INTERPRETATION 


The origin of the non-hydrostatic second zonal harmonic has been a source of speculation 
since the artificial satellite data determined a value of C 20 in excess of that predicted by 
hydrostatic theory. MacDonald (1963) postulated that AC 20 , the observed excess over 
the hydrostatic flattening due to the rotation of the earth, was a ‘fossil’ bulge remanent 
from the earth’s faster rotation rate 10 m yr ago. Maintenance of the fossil bulge required 
a deeper mantle viscosity of 10 26 poises, too high to allow mantle convection. Goldreich 
and Toomre (1969) showed that the ratio of the normalized non-hydrostatic AC 20 /C 22 
amounted to only 2.4 rather than the non-normalized ratio of 12 obtained by MacDonald 
and that this remaining excess can be explained as an artifact of the coordinate system. 

Wang (1966) suggested that the non-hydrostatic AC 20 resulted from Pleistocene glaciation; 
O’Connell (1971) showed that the non-hydrostatic bulge is 3 times too large to be accounted 
for by deglaciation; also, the results obtained above indicate that a field arising from excess 
polar anomalies will produce a high ratio of AC 20 /C 22 . A glacial origin of the third harmonic 
has also been discussed: Khan and O’Keefe (1974) compare the large negative Antarctic 
anomaly with the shrinkage of the ice cap in the Pliocene, although there is the difficulty of 
explaining the slow recovery from the melting of the ice load 4-5 m yr ago. This paper finds 
two additional difficulties: the mass deficiency of 2.3 X 10 22 g postulated in the above 
paper does not account more than for 1/3 of the observed zonal harmonic, and this sugges- 
tion cannot begin to explain the largest third harmonic coefficients: the GEM-9 field gives 
C 31 = 2.02 X 10~ 6 and S 33 = 1.41 X 10” 6 whereas theC 30 = .95 X 10~ 6 (Lerch, F. J. et 
al., 1977). 

Hide and Horai (1968) suggested the second harmonic could be attributed to undulations of 
the core-mantle boundary, but that the stresses required to support the undulations increased 
dramatically with the order of the harmonics taken. Seismic wave reflections from the 
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core-mantle interface have demonstrated that boundary undulations of more than 5-10 km 
do not exist (Engdahl and Johnson 5 1974;Buchbinder, 1968). 

Higbie and Stacey (1970) sought a single depth where all harmonic orders would produce 
comparable stresses. Their results were field-dependent, but they conclude that anomalies 
from depths below 900 km could not explain the fields. 

A widely employed method has considered the density anomalies as short wavelength fea- 
tures which are uncorrelated beyond brief distances. This method was outlined by Cook 
(1967) and employed by Guier and Newton (1965), Allan (1972), Khan (1977) and 
Lam beck (1976). Conclusions obtained from this technique were not always consistent: 
Allan favored a two-component field with effective depths of anomalous regions at 1700 km 
and 250 kin with a larger anomalous density at the lower layer, while Khan (1977) and 
Lambeck (1976) concluded that the measured anomalies arise from regions above 800-1000 
km. The mass point technique employed in this paper is theoretically equivalent to this 
technique and is used to evaluate the above conclusions. 

In the simulations described in the section above, it was found that considerable variations 
occur naturally in the low degree harmonics with changes of density models, changes either 
of type or of geographic position. Further, large sets of globally distributed random 
anomalies produce widely varying values of the second and third harmonics. Simulations of 
one mass point failed to locate a point on the globe which could simultaneously supply the 
larger third harmonic coefficients. Therefore it seems that these harmonics contain insuffi- 
cient information to establish any particular origin or site and that approach was abandoned 
in favor of analyses including all the harmonics and allowing the possibility that the density 
variations may be widely distributed. Implicit assumptions are that the density variations 
should have a reasonable magnitude and that radial distributions satisfying Kaula’s ‘rule of 
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thumb’ would be sought rather than special cases of geographic positions. For example, a 
model confining density anomalies to a thin surface layer (in the manner of the surface 
density representation) is possible but unlikely since a random set of surface anomalies 
would hot produce Kaula’s ‘rule of thumb’ 

The models considered in this paper fall into two broad classes 1) extended sources, that is 
long wavelength anomalies and 2) mass points, that is short wavelength anomalies. None of 
the models considered gravity anomalies originating from the core (below 3000 km) as it 
is difficult to suppose that a liquid outer core can support 'anoma lies of any substance. 

The impact of extended sources on the relative amplitude was calculated by Pollack (1973 
using an even density cap. Kaula (1977) computed the ‘degree variation’ (a fi ) using a Monte 
Carlo simulation of extended sources and determined that sources correlated within 20° 
best fit the observed field. He found that density anomalies must extend to 3000 km depth 
and that the values of decay length indicated that the anomalies must increase with depth. 
His model using the horizontal spacing y c ~ 20° best fit the exponent n of Kaula’s rule of 
thumb 1 O'" 5 /£ n but the closeness of the fit to the model significantly exceeded the observed 
fit. In particular, the mid-degree (6-11) spectra were substantially higher than predicted by 
Kaula’s rule of thumb. This can be understood by considering the convex curvature of the 
relative amplitude produced by extended sources, as shown by Pollack and demonstrated 
above, in contrast with the concavity of the relative amplitude of the observed field. Kaula 
also noted the excess scatter of the simulated models compared to the observed field. 

The results above indicate that extended sources are represented more heavily in the 2nd 
and 3rd harmonics (in comparison with mass points). Because the variation of the second 
and third harmonics from any model is large, it is difficult to assess the extent of the contri- 
bution of extended sources to these coefficients. However, these harmonics seem to have 
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larger values than would usually occur from models of random mass points so it is possible 
that extended sources contribute 1/3 or 1/4 of the magnitude of these coefficients. 

The abrupt decline in the relative amplitude of an extended source implies that such a feature 
will disappear when the long wavelengths are removed. Wagner (personal communication) 
has pointed out the Indian Ocean deficit is mainly represented in the low degree field (fi < 4). 
Gravity fields derived from local data types, rather than spherical harmonic representations 
may better measure extended sources. This possibility should be further studied as it is of 
interest in considering future techniques for gravity field recovery. Kaula (1963) noted that 
the observed field seemed to consist of incoherent short wavelength features; the above 
results demonstrate that the spherical harmonic representation of the gravity field measures 
such features in preference to extended features of equal mass. 

Simulations did not confirm that the field may be constructed of random short wavelength 
anomalies of equal density extending no deeper than 1000 km. This result is not in accord 
with earlier findings of Lambeck (1976) and Khan (1977) and thus further consideration of 
this discrepancy is warranted. Lambeck estimated that the density anomalies which con- 
tribute to the measured harmonics are located at less than 800 km depth, with the power in 
the harmonics above degree 6 attributed to heterogenieties below 300 or 400 km. In de- 
riving this result, he made two assumptions concerning the density variations; 1 ) that the 
power spectrum of the anomalies is the same for all layers of the spherical shells and 2) that 
the power spectrum may be represented by a white noise distribution. The first assumption 
is neither supported nor conclusively refuted by geophysical evidence. However, since there 
are sharp changes in density as a function of radius, it is quite possible that the average am- 
plitude of density anomalies may also vary with the radius. Thus the first assumption im- 
poses a serious constraint on the postulated mass distributions. 
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The second assumption, that the power spectrum of the density may be represented by a 
white noise spectrum, can be examined in greater detail. Lambeck expresses the anomalous 
density in spherical harmonics 

. if | ■ ■% f 

' = TT. < c im (P)P g ,„ w c° s + s cm (p)P Cm ( 6 ) sin mA) 

2=0 m=0 

and considers the power spectrum * 

e _ 

vf(A P )=2 (C L<p )+s Lw ) 

m=0 

where V^(Ap) is assumed to be a constant at all wavelengths. This assumed power may be 
compared with the spectral power of postulated density anomalies. The power of a mass 
point is proportional to 1/(2 £ + l) 1 ^, as discussed above. The spectrum of the point mass is 
a limiting case for non-organized elements: as Pollack (1973) stated, “a point mass displays 
the richest harmonic spectrum at all harmonic degrees.” That is, when the spectrum of a 
point mass is compared with extended mass distributions such as a spherical cap or spherical 
rectangle, that of the point mass declines less sharply with increasing degree. Further, the 
assumption of white noise implies that the amplitude of the shorter wavelengths are just as 
large as that of longer wavelengths, a geophysically unreasonable assumption. Although it is 
possible to construct a mass distribution approximating the power spectrum of a white 
noise spectrum, this is clearly a special event and may require unreasonable mass distribu- 
tions, If applied to a real power spectrum, the white noise assumption will exhibit a prefer- 
ence for a shallow depth of origin compared to a geophysically reasonable density 
distribution. 
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Application of this technique by Khan involved the minimizing of the residuals from a func- 
tion derived by Cook (1967) so that a fit to the slope could be obtained* This procedure 
reduces the sensitivity of the degree variance so that the result is less significant. It may be 
noted that all fields which are cumulative from the surface to an arbitrary depth D produce 
spherical harmonic coefficients which have the same order of magnitude as predicted by 
Kaula’s rule of thumb. The differences between fields become more discernible when the 
relative amplitude <I> C employed by Pollack (1973) is used. Further, applications of Cook’s 
technique usually involve the use of assumptions which predicate the answer in a complex 
manner. 

; l 

The results in this paper indicate that it is unlikely that the low harmonics arise from 
random mass anomalies near the surface because the amplitude does not decay as steeply as 
that of the observed field. Further, it seems difficult to match the observed amplitude with 
that arising from extended, sources as these produce an amplitude which as a convex slope 
rather than concaye, Thus, the observed gravitational field is not readily explained by 
assuming only shallow sources. 

Two models of anomalous density as a function of radius were found which did match the 
general characteristics of the observed relative amplitudes of the geopotential. One model is a 
modification of the ‘two-component’ model of Allan (1972), where two magnitudes of 
density anomalies in the upper and lower mantle are employed; the other model has density 
anomalies which continuously increase from the surface to the core-mantle boundary. These 
models represent extremes of a continuum of possible models which cannot be resolved in 
the external geopotential. Both models postulate deep mantle anomalies of larger magni- 
tude than the upper mantle and crustal anomalies. Since no models without this property 
were found to be capable of providing a satisfactory match to the observed relative amplitude, 
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it is concluded that deep mantle anomalies of average magnitude larger than those of the 
crust and upper mantle are a geophysical property of the earth. 

The modified two component model represents the lower mantle component, which forms 
the low degree harmonics, with layers containing an effective mass of 1 X 1CT 6 em each 100 
km from 1000 km to 2800 km depth (in contrast to a single layer at 1500 km in the Allan 
model). The upper component, which forms the high degree harmonics, is represented by a 
single layer of effective mass of 1 X 10" 6 em at 200 km depth. The lower component is 
larger by an order of magnitude since the upper component is the only mass in the upper 
1000 km whereas the lower component has mass layers every 100 km and by an additional 
factor of 4 to compensate for decrease in spherical area of a layer at 3000 km depth. Thus 
if the upper component is distributed throughout the 1000 km, this model requires an in- 
crease in anomalous density of a factor of 40 between the surface and the core; or, if the 
upper component is regarded as confined to the first 200 km, this model requires a region in 
which the anomalous density apparently vanishes. 

If the lower component is composed of 800 blocks of 5 X 5° in extent and 100 km in thick- 
ness, the effective mass of 1 X 10" 6 em at 1000 km depth may be converted to an average 
anomalous density <A p) = .01g/cm 3 , and at 2800 km (A p> = .04g/cm 3 . If, as this model 
may suggest, there is a region where little anomalous density exists, the geophysical implica- 
tion would be that the driving force of plate tectonics would be confined to motion in the 
upper mantle or crust. 

The modified two-component model has the additional deficiency of limiting the major 
contributions to the low harmonics to layers a few hundred km thick thus reducing the 
number N of density anomalies. Large N decreases the deviations of from the general 
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trend, which is desirable since the observed field demonstrates relatively little scatter in 
the <h c . j 

The continuously increasing anomalous density model is represented by mass points spaced 
100 km apart from the earth’s surface to the base of the mantle and increasing linearly from 
.25 X 10” 6 em at the surface to .95 X 10" 6 at 2800 km depth. At the surface, the average 
anomalous density <Ap> in a 5 X 5° block of J00 km thickness is found to be .0025g/cm 3 
and at 2800 km <Ap> = .03g/cm 3 , Thus a ratio of less than 1% in <Ap)/p can explain the 
gravitational field. 

This model also has the property that all layers contribute an equal amount to the second 
harmonic and nearly equal amounts to other low hannonics. Since the coefficients averaged 
over the largest number N of anomalies produce relative amplitudes with minimum scatter, 
this model will produce the smoothest field. 

The continuously increasing anomalous density model allows that sources of geophysical 
energy deep within the mantle (McKenzie, 1974) can be transmitted to the crust, providing 
an interior origin for plate tectonics (Davies, 1977). 

There is geophysical evidence for anomalous densities, including deep mantle anomalies. 
Seismological evidence for significant lateral variations in the lowest few hundred km of the 
mantle have been found by Julian and Sengupta (1973), and Jordan and Lynn (1974). 
Recently, a global inversion of travel times of first arrivals of compressional waves has been 
obtained by Dziewonski et al. (1977). The earth was divided into 5 spherical layers of 
depths from 0-670 km, 670-1 100 km, 1100-1500 km, and 2200-2800 km. Each layer was 
divided into 30 blocks with dimensions of 60° longitude and 36° latitude. The data base 
consisted of 700,000 travel time residuals formed by differencing the observed travel time 
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with the travel time computed from tables such as the Jeffrey-Bullen travel time tables. An 
inversion procedure was used to recover the anomalous velocity within each block from the 
travel time residual. Two discretization patterns for blocks were used to obtain two solu- 
tions for the seismic velocity deviations within a block; the two solutions for the three lower 
layers were similar although the two layers above 1 1 00 km showed little correlation between 
the solutions. These solutions were converted in spherical harmonic expressions for the seis- 
mic velocity deviations taken through degree and order 3. The coefficients of the seismic 
velocity deviation of the spherical harmonic expression were found to be correlated signifi- 
cantly with the ten largest significant coefficients of the geopotential. In fact, the authors 
found that the gravity and seismic fields could be matched exactly within the limits of the 
accuracy of the seismic velocity deviation data. Contour maps derived from the spherical 
harmonic model of the velocity anomalies indicate sizable features exist in all layers of the 
mantle with the largest deviations at the base of the mantle. 

Maximum seismic deviations in a 60° block in the deepest region amount to 60 m/s. If the 
ratio Au/Ap = -6(km/s)/(g/cm 3 ) derived from comparing seismic and gravity field coeffi- 
cients by Dziewonski et al. is employed, then Ap = .01 g. This is consistent with the above 
values of Ap, derived from consideration of gravity alone. 

The existence of deep-seated lateral seismic inhomogeneities provides independent support 
for deep density anomalies, in addition to the compelling evidence obtained from the 
simulations above. 
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IV. CONCLUSIONS 


Two models of radial distribution of anomalous density were found that can replicate the 
energy spectrum of the observed gravitational field. Both models contained deep anomalous 
densities of larger magnitude than shallow anomalies. These deep anomalies extend to the 
base of the mantle in both models. No model of anomalous density as a function of radius 
can replicate the properties of the observed field as a function of degree if the density anom- 
alies exist only above 1000 km depth. Therefore it is probable that deep and large density 
anomalies form a significant portion of the earth’s geopotential contained in the low degree and 
order spherical harmonic coefficients. The minimum depth at which there must be large anom- 
alies is 1500 km, but very large anomalies would be required unless the anomalies are extended 
deeper. Anomalies extending to the base of the mantle with average anomalous densities of 
<Ap> ~ .03g/cm 3 can represent the relative amplitude as a function of degree. Lateral varia- 
tions in seismic travel times support the concept of a heterogeneous lower mantle. This con- 
cept is also in accord with the notable lack of correlation of low degree harmonics with sur- 
face features (Kaula, 1967). 

The external geopotential expressed in spherical harmonics preferentially measures uncorre- 
lated short wavelength features over long wavelength features, with the exception of the very 
low degrees where laterally extended features may be significant. The summation of several 
short wavelength features produces a noticeable scatter of the energy spectrum about the 
trend line produced by one similar feature. Several thousand short wavelength features may 
be necessary to produce a simulated field with a degree of scatter in the energy spectrum 
comparable to the relatively small scatter observed. 

The observed value of the energy of the second harmonic is compatible in magnitude with . 
the models, although there is an excess of ~10% compared to most of the models which 
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well match the degrees above two. This possible excess can be attributed either to scatter, 
which is largest in the second degree, or to the contributions of extended features. It is not 
possible to recover any significant geophysical information from the non-hydrostatic second 
zonal harmonic, per se. 

The two models of radial anomalous density able to represent the decline of the energy 
spectrum are extremes of a continuum of satisfactory models. The monotonically increas- 
ing anomalous density model postulates a continuously increasing density, where the effec- 
tive anomalous mass/spherical layer has increased by a factor of 4 from the surface to the 
base of the mantle to about 1 X 10" 6 earth masses at 2800 km depth. The modified two- 
component model postulates an effective mass of 1 X 10" 6 em in the region between 0 and 
1000 km depth and an effective mass of 1 X 10~ 6 em on spherical layers placed every 100 
km between 1000 and 2800 km depth. Discontinuous changes in anomalous density where 
the radial density is also discontinuous are reasonable, but if so, such changes cannot be re- 
solved in the gravity observations. 

The monotonically increasing anomalous density model may be related to recent discussions 
of the geophysics of the mantle. In a review, O’Connell (1977) showed that the data indi- 
cates there is no compelling reason for mantle convection to be confined to the upper 
mantle and there is evidence that mantle convection extends to the lower mantle. Recent 
studies have found a uniform viscosity for the mantle (Dicke, 1969; O’Connell, 1971; 
Cathles, 1975). A review of thermal convection calculations together with other evidence 
concludes that unless there are large viscosity differences (~10 4 ) the flow induced by plate 
tectonic motion will penetrate into the lower mantle (Davies, 1977). Since the results above 
demonstrate that the earth’s gravity field is best explained by the existence of strong deep 
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mantle density anomalies, it seems reasonable to suppose that lower mantle convection is 
the agent of plate tectonic motion rather than the consequence. 

The anomalies in the crust and upper mantle predominate in the higher degree field. 

With the advent of high degree and order gravity fields, it became possible to correlate sur- 
face features with gravity anomalies. Kaula (1969, 1972) found strong associations of 
anomalies with extensive Quaternary volcanism and orogeny, active ocean rise and ocean 
basins. Other correlations with near surface features include long-wavelength elevations 
associated with cooling plates spreading from ocean rises (Sclater et al., 1975) and Pacific 
ocean anomalies associated with aesthenospheric convection rolls (Marsh and Marsh, 1976; 
Richter and Parsons, 1975). These features predominate in the higher degrees; however, as 
Rapp (1977) showed, the degree at which an extended feature is energetic is not related in a 
simple manner to the angular size of the feature. 

In the monotonically increasing model, every 100 km thick shell between the surface and 
the core of the earth provides a contribution of ~3% to the energy of the non-hydrostatic 
second degree harmonic. In this event, the origins of the harmonic coefficients are every- 
where and gravitational data must be combined with other geophysical data in order to 
interpret lateral density heterogeneities. 
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